function [perm_quadrant_summary_table, perm_cubic_summary_table] = generate_cross_validation_permutation_100_table(pm,split,reweight,covariates,cost,SMC,s_winter,target_wave,sample_wave,fbase,savedir_figures)
% This function summarizes out-of-sample performance of rules in the pooled
% sample.

% Inputs: 
% (1) pm: number permutations
% (2) split: share of data in the testing set
% (3) reweight: use density ratio to reweight waves or not
% (4) covariates: covariate specification
% (5) cost: private cost savings, social cost savings, or energy savings
% (6) SMC: =1 use social marginal cost; =0 use retail electricity price
% (7) s_winter: ='' since include all baseline periods
% (8) target_wave: =0 in pooled sample
% (9) sample_wave: =0 in pooled sample
% (10) fbase: string used to indicate the propensity score and baseline
% months specification for output table 
% (11) savedir_figures: directory to save the figures

% Outputs: 
% (1) perm_quadrant_summary_table: table summarizes out-of-sample savings
% for quadrant rules
% (2) perm_cubic_summary_table: table summarizes out-of-sample savings for 
% % cubic rules  

%% Summarize quadrant ewm savings and percentage of treatment
ewm_training_quadrant_perm=nan(pm,1);
ewm_testing_quadrant_perm=nan(pm,1);
ewm_training_quadrant=cell(size(cost,2),size(covariates,2));
ewm_testing_quadrant=cell(size(cost,2),size(covariates,2));
percent_training_quadrant_perm=nan(pm,1);
percent_testing_quadrant_perm=nan(pm,1);
percent_training_quadrant=cell(size(cost,2),size(covariates,2));
percent_testing_quadrant=cell(size(cost,2),size(covariates,2));
att_testing_quadrant_perm=nan(pm,1);
att_training_quadrant_perm=nan(pm,1);
att_testing_quadrant=cell(size(cost,2),size(covariates,2));
att_training_quadrant=cell(size(cost,2),size(covariates,2));
matrix_quadrant = cell(size(cost,2),size(covariates,2));

s_wave=sprintf('%d%d',sample_wave,target_wave);

for j=1:size(covariates,2) % covariate
    for i=1:size(cost,2) % cost
        if (cost{i})
            if (SMC{i})
                s_cost = 'SMC';
            else
                s_cost = 'PMC';
            end
        else
            s_cost='kwh';
        end
            matrix_quadrant{i,j} = strcat(covariates{j},'_',s_cost); % save quadrant specification matrix

            for p=1:pm % number repeated runs
                s_fold=sprintf('%d',p);
                s_split=sprintf('split%.1f',split);

                if reweight==0
                    filename_coefs=sprintf('coef_cv_quadrant_%s_%s%s_%s_%s_perm%s.mat',covariates{j},s_cost,s_winter,s_wave,s_split,s_fold);
                else
                    filename_coefs=sprintf('coef_cv_quadrant_%s_%s%s_%s_%s_perm%s_reweight.mat',covariates{j},s_cost,s_winter,s_wave,s_split,s_fold);
                end

                filename = sprintf('%s_%s',fbase,filename_coefs);
                S=load(filename);
                ewm_testing_quadrant_perm(p,1)=sum(S.gu_target_set.*S.in_Ghat)/S.n_target_set;
                ewm_training_quadrant_perm(p,1)=sum(S.gu_w.*S.in_Ghat)/S.n_sample_set;
                percent_testing_quadrant_perm(p,1)=sum(S.nw_target_set.*S.in_Ghat)/S.n_target_set;
                percent_training_quadrant_perm(p,1)=sum(S.nw_sample_set.*S.in_Ghat)/S.n_sample_set;
                att_testing_quadrant_perm(p,1)=S.scaled_att_target_set;
                att_training_quadrant_perm(p,1)=S.scaled_att_sample_set;

            end
        ewm_training_quadrant{i,j}=ewm_training_quadrant_perm;    
        ewm_testing_quadrant{i,j}=ewm_testing_quadrant_perm;
        percent_training_quadrant{i,j}=percent_training_quadrant_perm;
        percent_testing_quadrant{i,j}=percent_testing_quadrant_perm;
        att_testing_quadrant{i,j}=att_testing_quadrant_perm;
        att_training_quadrant{i,j}=att_training_quadrant_perm;
    end
end  

%% Summarize univariate ewm savings and percentage of treatment
ewm_training_onedim_perm=nan(pm,1);
ewm_testing_onedim_perm=nan(pm,1);
ewm_training_onedim=cell(size(cost,2),size(covariates,2));
ewm_testing_onedim=cell(size(cost,2),size(covariates,2));
percent_training_onedim_perm=nan(pm,1);
percent_testing_onedim_perm=nan(pm,1);
percent_training_onedim=cell(size(cost,2),size(covariates,2));
percent_testing_onedim=cell(size(cost,2),size(covariates,2));
att_testing_onedim_perm=nan(pm,1);
att_training_onedim_perm=nan(pm,1);
att_testing_onedim=cell(size(cost,2),size(covariates,2));
att_training_onedim=cell(size(cost,2),size(covariates,2));
matrix_onedim = cell(size(cost,2),1);

s_wave=sprintf('%d%d',sample_wave,target_wave);

for j=1:size(covariates,2) % covariate
    for i=1:size(cost,2) % cost
        if (cost{i})
            if (SMC{i})
                s_cost = 'SMC';
            else
                s_cost = 'PMC';
            end
        else
            s_cost='kwh';
        end
        matrix_onedim{i,1} = strcat(s_cost); % save univariate specification matrix 

            for p=1:pm 
                s_fold=sprintf('%d',p);
                s_split=sprintf('split%.1f',split);

                if reweight==0
                    filename_coefs=sprintf('coef_cv_onedim_%s_%s%s_%s_%s_perm%s.mat',covariates{j},s_cost,s_winter,s_wave,s_split,s_fold);
                else
                    filename_coefs=sprintf('coef_cv_onedim_%s_%s%s_%s_%s_perm%s_reweight.mat',covariates{j},s_cost,s_winter,s_wave,s_split,s_fold);
                end

                filename = sprintf('%s_%s',fbase,filename_coefs);
                S=load(filename);
                ewm_testing_onedim_perm(p,1)=sum(S.gu_target_set.*S.in_Ghat)/S.n_target_set;
                ewm_training_onedim_perm(p,1)=sum(S.gu_w.*S.in_Ghat)/S.n_sample_set;
                percent_testing_onedim_perm(p,1)=sum(S.nw_target_set.*S.in_Ghat)/S.n_target_set;
                percent_training_onedim_perm(p,1)=sum(S.nw_sample_set.*S.in_Ghat)/S.n_sample_set;
                att_testing_onedim_perm(p,1)=S.scaled_att_target_set;
                att_training_onedim_perm(p,1)=S.scaled_att_sample_set;
            end
        ewm_training_onedim{i,j}=ewm_training_onedim_perm;    
        ewm_testing_onedim{i,j}=ewm_testing_onedim_perm;
        percent_training_onedim{i,j}=percent_training_onedim_perm;
        percent_testing_onedim{i,j}=percent_testing_onedim_perm;
        att_testing_onedim{i,j}=att_testing_onedim_perm;
        att_training_onedim{i,j}=att_training_onedim_perm;
    end
end  

%% Summarize cubic ewm savings and percentage of treatment
ewm_training_cubic_perm=nan(pm,1);
ewm_testing_cubic_perm=nan(pm,1);
ewm_training_cubic=cell(size(cost,2),size(covariates,2));
ewm_testing_cubic=cell(size(cost,2),size(covariates,2));
percent_training_cubic_perm=nan(pm,1);
percent_testing_cubic_perm=nan(pm,1);
percent_training_cubic=cell(size(cost,2),size(covariates,2));
percent_testing_cubic=cell(size(cost,2),size(covariates,2));
att_testing_cubic_perm=nan(pm,1);
att_training_cubic_perm=nan(pm,1);
att_testing_cubic=cell(size(cost,2),size(covariates,2));
att_training_cubic=cell(size(cost,2),size(covariates,2));


for j=1:size(covariates,2) % covariate
    for i=1:size(cost,2) % cost
        if (cost{i})
            if (SMC{i})
                s_cost = 'SMC';
            else
                s_cost = 'PMC';
            end
        else
            s_cost='kwh';
        end

            for p=1:pm 
                s_fold=sprintf('%d',p);
                s_split=sprintf('split%.1f',split);
                
                if reweight==0
                    filename_coefs=sprintf('coef_cv_cubic_%s_%s%s_%s_%s_perm%s.mat',covariates{j},s_cost,s_winter,s_wave,s_split,s_fold);
                else
                    filename_coefs=sprintf('coef_cv_cubic_%s_%s%s_%s_fold%s_reweight.mat',covariates{j},s_cost,s_winter,s_wave,s_fold);
                end

                filename = sprintf('%s_%s',fbase,filename_coefs);
                S=load(filename);
                ewm_testing_cubic_perm(p,1)=nanmean(S.g_target_set.*S.in_Ghat)*S.Yscale_target_set;
                in_Ghat_sample=(S.X_sample_set*S.beta>0);
                ewm_training_cubic_perm(p,1)=nanmean(S.g_sample_set.*in_Ghat_sample)*S.Yscale_sample_set;
                percent_testing_cubic_perm(p,1)=mean(S.in_Ghat);
                percent_training_cubic_perm(p,1)=mean(in_Ghat_sample);
                att_testing_cubic_perm(p,1)=S.scaled_att_target_set;
                att_training_cubic_perm(p,1)=S.scaled_att_sample_set;
            end
        ewm_training_cubic{i,j}=ewm_training_cubic_perm;    
        ewm_testing_cubic{i,j}=ewm_testing_cubic_perm;
        percent_training_cubic{i,j}=percent_training_cubic_perm;
        percent_testing_cubic{i,j}=percent_testing_cubic_perm;
        att_testing_cubic{i,j}=att_testing_cubic_perm;
        att_training_cubic{i,j}=att_training_cubic_perm;
    end
end  

%% organize treatment share and savings into a table - quadrant & univariate
perm_quadrant_summary_table=nan(5,5);
perm_quadrant_summary_table=array2table(perm_quadrant_summary_table,'VariableNames',{'Rule','Covariates','kwh','PMC','SMC'});
format short g
perm_quadrant_summary_table.(1) = {'RCT';'EWM-univariate';'EWM-quadrant';'EWM-quadrant';'EWM-quadrant'};
perm_quadrant_summary_table.(2) = {'\';'\';covariates{1};covariates{2};covariates{3}}; % covariate column
perm_quadrant_summary_table.(3) = {'\';'\';'\';'\';'\'}; 
perm_quadrant_summary_table.(4) = {'\';'\';'\';'\';'\'}; 
perm_quadrant_summary_table.(5) = {'\';'\';'\';'\';'\'}; 

for j=1:size(covariates,2) % covariate
    covariate = covariates{j};
    for i=1:size(cost,2) % cost
        if (cost{i})
            if (SMC{i})
                s_cost = 'SMC';
            else
                s_cost = 'PMC';
            end
        else
            s_cost='kwh';
        end
         % fill in quadrant savings 
        table_column_idx=find(strcmp(perm_quadrant_summary_table.Properties.VariableNames',s_cost)); % find the column index in table that corresponds to s_cost
        table_row_idx=find(strcmp(covariate,perm_quadrant_summary_table.Covariates)); % find the row index in table that corresponds to covariate 
        find_in_matrix = strcmp(strcat(covariate,'_',s_cost),matrix_quadrant); % find corresponding number in the savings or percent matrix, according to covariate_scost 
        perm_quadrant_summary_table{table_row_idx,table_column_idx} = {sprintf('%0.3f\x00B1%0.2f',mean(ewm_testing_quadrant{find_in_matrix}),var(ewm_testing_quadrant{find_in_matrix})^(1/2))}; % extract numbers according to index match
        
        % fill in univariate savings 
        table_column_idx=find(strcmp(perm_quadrant_summary_table.Properties.VariableNames',s_cost));
        table_row_idx=find(strcmp('EWM-univariate',perm_quadrant_summary_table.Rule));
        find_in_matrix = strcmp(strcat(s_cost),matrix_onedim);
        perm_quadrant_summary_table{table_row_idx,table_column_idx} = {sprintf('%0.3f\x00B1%0.2f',mean(ewm_testing_onedim{find_in_matrix}),var(ewm_testing_onedim{find_in_matrix})^(1/2))};
        
        % fill in RCT savings 
        table_column_idx=find(strcmp(perm_quadrant_summary_table.Properties.VariableNames',s_cost));
        table_row_idx= 1; % RCT stays in the top row
        find_in_matrix = strcmp(strcat(s_cost),matrix_onedim);
        perm_quadrant_summary_table(table_row_idx,table_column_idx) = {sprintf('%0.3f\x00B1%0.2f',mean(att_testing_onedim{find_in_matrix}),var(att_testing_onedim{find_in_matrix})^(1/2))};
    end
end

%% organize treatment share and savings into a table - cubic

perm_cubic_summary_table=nan(4,5);
perm_cubic_summary_table=array2table(perm_cubic_summary_table,'VariableNames',{'Rule','Covariates','kwh','PMC','SMC'});
format short g
perm_cubic_summary_table.(1) = {'RCT';'EWM-cubic';'EWM-cubic';'EWM-cubic'};
perm_cubic_summary_table.(2) = {'\';covariates{1};covariates{2};covariates{3}}; % covariate column
perm_cubic_summary_table.(3) = {'\';'\';'\';'\'}; 
perm_cubic_summary_table.(4) = {'\';'\';'\';'\'}; 
perm_cubic_summary_table.(5) = {'\';'\';'\';'\'}; 

for j=1:size(covariates,2) % covariate
    covariate = covariates{j};
    for i=1:size(cost,2) % cost
        if (cost{i})
            if (SMC{i})
                s_cost = 'SMC';
            else
                s_cost = 'PMC';
            end
        else
            s_cost='kwh';
        end
        
        % fill in cubic savings 
        table_column_idx=find(strcmp(perm_cubic_summary_table.Properties.VariableNames',s_cost)); % find the column index in table that corresponds to s_cost
        table_row_idx=find(strcmp(covariate,perm_cubic_summary_table.Covariates)); % find the row index in table that corresponds to covariate 
        find_in_matrix = strcmp(strcat(covariate,'_',s_cost),matrix_quadrant); % find corresponding number in the savings or percent matrix, according to covariate_scost 
        perm_cubic_summary_table{table_row_idx,table_column_idx} = {sprintf('%0.3f\x00B1%0.2f',mean(ewm_testing_cubic{find_in_matrix},'omitnan'),var(ewm_testing_cubic{find_in_matrix},'omitnan')^(1/2))}; % extract numbers according to index match
    
        % fill in RCT savings 
        table_column_idx=find(strcmp(perm_cubic_summary_table.Properties.VariableNames',s_cost));
        table_row_idx= 1; % RCT stays in the top row
        find_in_matrix = strcmp(strcat(s_cost),matrix_onedim);
        perm_cubic_summary_table(table_row_idx,table_column_idx) = {sprintf('%0.3f\x00B1%0.2f',mean(att_testing_onedim{find_in_matrix},'omitnan'),var(att_testing_onedim{find_in_matrix},'omitnan')^(1/2))};
    end
end


end